def make_ms_frame(msname:str, ant_config:str, spw_params=None, verbose=False):
"""
Construct an empty MeasurementSet that has the desired observation setup.
Args:
msname: MS to create
ant_config: a telescope configuration from casadata (alma/simmos) (for
example alma/simmos/aca.cycle8.cfg)
swp_params: parameters such as SPW name, number of channels, frequency
"""
## Open the simulatorFunctions
sm.open(ms=msname);
## Read/create an antenna configuration.
## Canned antenna config text files are located here: <casadata>/alma/simmos/*cfg
antennalist = os.path.join(ctsys.resolve("alma/simmos"), ant_config)
if verbose:
print(f'Using antenna list file: {antennalist}')
## Fictitious telescopes can be simulated by specifying x, y, z, d, an, telname, antpos.
## x,y,z are locations in meters in ITRF (Earth centered) coordinates.
## d, an are lists of antenna diameter and name.
## telname and obspos are the name and coordinates of the observatory.
mysu = simutil.simutil()
(x, y, z, d, an, an2, telname, obspos) = mysu.readantenna(antennalist)
## Set the antenna configuration
metool = measures()
sm.setconfig(
telescopename=telname,
x=x,
y=y,
z=z,
dishdiameter=d,
mount=['alt-az'],
antname=an,
coordsystem='local',
referencelocation=metool.observatory(telname)
)
## Set the polarization mode (this goes to the FEED subtable)
sm.setfeed(mode='perfect X Y', pol=['']); # X Y / R L
## Set the spectral window and polarization (one data-description-id).
## Call multiple times with different names for multiple SPWs or pol setups.
sm.setspwindow(
spwname=spw_params['name'],
freq=spw_params['freq'],
deltafreq='0.1GHz',
freqresolution='0.2GHz',
nchannels=spw_params['nchan'],
stokes='XX YY'
)
## Setup source/field information (i.e. where the observation phase center is)
## Call multiple times for different pointings or source locations.
source_name = 'simulated_source'
sm.setfield(
sourcename=source_name,
sourcedirection=metool.direction(
rf='J2000',
v0='19h59m28.5s',
v1='+40d44m01.5s'
)
)
## Set shadow/elevation limits (if you care). These set flags.
sm.setlimits(shadowlimit=0.01, elevationlimit='1deg');
## Leave autocorrelations out of the MS.
sm.setauto(autocorrwt=0.0);
## Set the integration time, and the convention to use for timerange specification
## Note : It is convenient to pick the hourangle mode as all times specified in sm.observe()
## will be relative to when the source transits.
sm.settimes(
integrationtime='2s',
usehourangle=True,
referencetime=metool.epoch('UTC','2021/10/14/00:01:02')
)
## Construct MS metadata and UVW values for one scan and did
## Call multiple times for multiple scans.
## Call this with different sourcenames (fields) and spw/pol settings as defined above.
## Timesteps will be defined in intervals of 'integrationtime', between starttime and stoptime.
sm.observe(
sourcename=source_name,
spwname=spw_params['name'],
starttime='0',
stoptime='+8s'
)
## Close the simulator
sm.close()
## Unflag everything (unless you care about elevation/shadow flags)
flagdata(vis=msname, mode='unflag', flagbackup=False)
def add_gaussian_noise(msname:str, noise_level_sigma='0.1Jy'):
"""
Adds Gaussian random noise using simulator / corrupt
Args:
ms_name: MS to modify
noise_level_sigma: noise sigma as used in the simulator corrupt function
"""
try:
sm.openfromms(msname)
sm.setseed(4321)
sm.setnoise(mode='simplenoise', simplenoise=noise_level_sigma)
sm.corrupt()
finally:
sm.close()
def add_point_source_component(msname, cl_name = 'sim_point_source.cl', freq=None,
flux=5.0, spectrumtype='constant', index=-1):
"""
Adds a point source to the MS
Args:
ms_name: MS to modify
freq:
spectrumtype:
flux: In Jy, as used in componentlist.addcomponent
"""
make_point_source_comp_list(
cl_name=cl_name,
freq=freq,
flux=flux,
spectrumtype=spectrumtype,
index=index
)
sim_from_comp_list(msname, cl_name)
shutil.rmtree(cl_name)
def make_point_source_comp_list(cl_name:str, freq:str, flux:str, spectrumtype:str, index:int,
direction='J2000 19h59m28.5s +40d44m01.5s',
fluxunit='Jy', shape='point', label='sim_point_source'):
"""
Makes a component list file with a point source
Args:
cl_name: name of the component list file to create
freq: freq quantity (with units) as used in the component list tool
flux: flux, units are assumed in Jy
spectrumtype: as used in component list tool: constant/spectral index
index: spectral index
"""
try:
cl.addcomponent(
dir=direction,
flux=flux,
fluxunit=fluxunit,
freq=freq,
shape=shape, # shape='gaussian',
# majoraxis="5.0arcmin", minoraxis='2.0arcmin',
# polarization='linear', / 'Stokes'
# spectrumtype:'spectral index' / 'constant'
spectrumtype=spectrumtype,
index=index,
label=label)
cl.rename(filename=cl_name)
finally:
cl.close()
def sim_from_comp_list(msname:str, cl_name:str):
"""
Updates the MS visibilities using simulator.predict to add
components from the components list
Args:
ms_name: MS to modify
cl_name: name of components list file to simulate
"""
try:
sm.openfromms(msname)
sm.predict(complist=cl_name, incremental=False)
finally:
sm.close()
def add_spectral_line(msname:str, line:'numpy.array', chan_range=[60, 86], amp_factor=None):
"""
Adds a spectral line as a Gaussian function in the range of channels given
Args:
ms_name: MS to modify
line: function to produce spectral line. ex: fn(x, mu, sigma)
chan_range: list with indices of the first and last channel
amp_factor: factor to multiply the peak height / flux density
"""
try:
tbtool = table()
tbtool.open(msname, nomodify=False)
data = tbtool.getcol('DATA')
if not amp_factor:
amp_factor = 1.0/np.max(line)
data[:,chan_range[0]:chan_range[1],:] += amp_factor * (1+0j) * line.reshape((1, len(line), 1))
tbtool.putcol('DATA', data)
finally:
tbtool.done()
def add_polynomial_continuum(msname:str, pol_coeffs:list, nchan:int, amp_factor=1.0):
"""
Update MS visibilities adding a polynomial evaluated for all channels.
Args:
ms_name: MS to modify
pol_coeff: polynomial coefficients, evaluated [0.5, 0.25] => 0.5x + 0.25
nchan: number of channels in the SPW (x axis to eval polynomial)
"""
try:
tbtool = table()
tbtool.open(msname, nomodify=False)
data = tbtool.getcol('DATA')
x = np.linspace(0, 1, nchan)
polynomial = np.polyval(pol_coeffs, x)
# Add same polynomial to real and imag part
data += amp_factor * (1+1j) * polynomial.reshape((1, len(polynomial),1))
tbtool.putcol('DATA', data)
finally:
tbtool.done()
